Julio Ramírez Guerrero, Julián Román Camacho, Silvestre Ruano Rodríguez, Manuel Racero de la Rosa, Rafael Rubio Ramos.
Ante la alarmante problemática climática actual, se pretende encontrar plantas que resistan condiciones de estrés abiótico. Los flavonoides (metabolito secundario sintetizados por la planta cuando se encuentra sometida a estrés) pueden ser de gran ayuda para determinar la relación planta-microbioma , y ver como esta mejora la resistencia al estrés hídrico. Esta relación está bien estudiada en especies de bacterias nodulantes, pero no tanto en bacterias no nodulantes como las de la familia Aeromonadaceae. Por tanto, la motivación del estudio es comprobar el comportamiento de plantas de Arabidopsis thaliana en presencia/ausencia de Aeromonas sp (He et al., 2022).
Comprobar la expresión diferencial en Arabidopsis thaliana, en presencia/ausencia de Aeromonas sp. y condiciones de sequía.
Para nuestro diseño experimental hemos usado como referencia el número de acceso GSE184872 de la web GEO (Gene Expression Omnibus).
Tenemos una serie de puntos que describen las condiciones de nuestro estudio, según el tratamiento que hemos aplicado a Arabidopsis thaliana, planta que es el organismo modelo de este estudio.
Sin tratamientos (Controles):
Tratamiento en sequía:
Tratamiento con Aeromonas:
Tratamientos en sequía y con Aeromonas:
A continuación, se muestra la serie de pasos que se han seguido para realizar el análisis matemático/computacional de los datos RNA-seq.
En primer lugar, se prepara el espacio de trabajo en linux, aquí se toman los datos del artículo de referencia y se descargan mediante sbatch. Esto proporcionará los archivos tanto .gz como .gtf. Mediante un fastq-dump se convierten los archivos SRA en FASTQ, posteriormente se realizará el análisis de calidad con fastqc de los fastq.gz y se revisará que esten correctamente analizados sin contaminaciones u otra alteración.
Seguido, se ejecutará hisat2 en los fastq.gz y así, se mapearán las lecturas de secuencias cortas contra el genoma de referencia, dando archivos .sam que pesan demasiado y son ineficientes a la hora de manipular los datos, por lo que se pasan a .bam mediante samtools obteniendo así archivos de menor peso. También se generarán los índices. Con los alineamientos de hisat2 se realiza un ensamblado de los transcritos de nuestro .bam con la herramienta stringtie. Una vez ensambladas cada una de las muestras en los correspondientes .gtf es necesario recopilar el transcriptoma completo, por ello se unen los .gtf con merge stringtie. Y una vez se obtiene el transcriptoma completo se cuantifica la expresión génica.
Se continua el análisis en el programa R. Con la librería ballgown se leen los datos de expresión obtenidos anteriormente, se establecen las 4 condiciones del análisis, cada una con sus 3 réplicas. Se extraen los niveles de expresión tanto por FPKM y TPM, luego, se normalizan los datos y se comparan los boxplot antes y después de normalizar. Se comprueba la diferencia entre réplicas con un scatterplot. Con el análisis de PCA y clustering jerárquico se verifican las diferentes distribuciones de los grupos y su variabilidad. Con la matriz de expresión se agrupan los datos y se comprueban su distribución con scatterplots que representan las diferentes condiciones y volcano plots para estudiar la magnitud biológica de las diferencias en la expresión génica entre diferentes condiciones experimentales. Y por último, un enriquecimiento funcional para destacar que procesos biológicos están más representados.
En el gráfico correspondiente al análisis de calidad de la muestra 1, los resultados obtenidos son de alta calidad, ya que los boxplot del gráfico se mantienen en todo momento en el área coloreada en verde. En el resto de muestras los resultados son similares, manteniéndose la alta calidad de las mismas.
Es necesario realizar un pre-procesamiento de los datos para poder comparar los datos de calidad de distintas muestras. Se empleó un algoritmo de normalización para normalizar los datos y poder realizar el análisis adecuadamente.
library(ballgown)
pheno.data <- read.csv("pheno_data.csv")
pheno.data
## sample aeromonas drought
## 1 sample01 ctl ctl
## 2 sample02 ctl ctl
## 3 sample03 ctl ctl
## 4 sample04 ctl drought
## 5 sample05 ctl drought
## 6 sample06 ctl drought
## 7 sample07 h1 ctl
## 8 sample08 h1 ctl
## 9 sample09 h1 ctl
## 10 sample10 h1 drought
## 11 sample11 h1 drought
## 12 sample12 h1 drought
bg.data <- ballgown(dataDir = ".", samplePattern = "sample", pData=pheno.data)
bg.data
## ballgown instance with 54013 transcripts and 12 samples
sampleNames(bg.data)
## [1] "sample01" "sample02" "sample03" "sample04" "sample05" "sample06"
## [7] "sample07" "sample08" "sample09" "sample10" "sample11" "sample12"
gene.expression <- gexpr(bg.data)
dim(gene.expression)
## [1] 32833 12
gene.names <- rownames(gene.expression)
colnames(gene.expression) <- c("ctl_water1","ctl_water2","ctl_water3","ctl_drought1","ctl_drought2","ctl_drought3", "h1_water1","h1_water2","h1_water3","h1_drought1","h1_drought2","h1_drought3")
gene.expression.1 <- gene.expression + 1
write.table(x = gene.expression.1,file = "tarea1_gene_expression.tsv",
quote = F,row.names = F,
sep = "\t")
boxplot(gene.expression, outline=F,col=rainbow(12),ylab="Gene Expression (FPKM)",
cex.lab=1.5)
boxplot(log2(gene.expression.1), outline=F,col=rainbow(12),
ylab="log2(Gene Expression)",
cex.lab=1.5)
Se hace el análisis de la expresión de los datos por TPM, para comprobar que los resultados son iguales que los realizados por FPKM:
tpm <- gene.expression.1 / colSums(gene.expression.1) * 1e6
## Guardar los datos de expresión génica en formato TPM
write.table(x = tpm, file = "tarea1_gene_expression_tpm.tsv", quote = F, row.names = F, sep = "\t")
## Representación de la distribución global de la expresión génica en TPM:
boxplot(tpm, outline = F, col = rainbow(12), ylab = "Gene Expression (TPM)", cex.lab = 1.5)
La normalización de datos es necesaria ya que durante la preparación de muestras pueden producirse errores que alteren los datos reales que se deberían de obtener. De este modo, corregimos el sesgo producido durante la preparación de muestras y facilita la comparación de resultados entre distintas muestras. Además, la normalización de datos facilita la detección de diferencias significativas en la expresión génica, y previene interpretaciones erróneas porque la normalización elimina las diferencias aparentes que podrían ser resultado de variaciones técnicas en lugar de biológicas
library(NormalyzerDE)
design <- data.frame(sample=colnames(gene.expression),
group=c(rep("ctl_water",3),rep("ctl_drought",3),rep("h1_water",3),rep("h1_drought",3)))
write.table(x = design,file = "normalyzer_design.tsv",quote = F,row.names = F,
sep = "\t")
normalyzer(jobName = "Tarea1",designPath = "normalyzer_design.tsv",
dataPath = "tarea1_gene_expression.tsv",outputDir = ".")
normalized.gene.expression <- read.table(file="Tarea1/Quantile-normalized.txt", header=T)
rownames(normalized.gene.expression) <- gene.names
Según esta evaluación VSN es la mejor técnica de normalización porque es la que más reduciría el ruido.
Sin embargo, según el coeficiente de correlación de Pearson VSN no es una técnica adecuada. Por lo que tomamos Quantil como técnica de normalización.
boxplot(normalized.gene.expression, outline=F,col=rainbow(12),
ylab="log2(FPKM + 1)",cex.lab=1.5)
Los diagramas salen alineados, indicando que están en la misma escala y son comparables, por lo tanto, se confirma que se ha realizado la normalización de datos.
Scatterplots para comparar réplicas
par(mfrow=c(3,4), mai = c(0.4, 0.4, 0.4, 0.4), mgp=c(2,1,0))
plot(x = normalized.gene.expression[,"ctl_water1"],
y = normalized.gene.expression[,"ctl_water2"],
pch=19,col="black",xlab="ctl_water1",ylab="ctl_water2",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"ctl_water1"],
normalized.gene.expression[,"ctl_water2"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"ctl_water1"],
y = normalized.gene.expression[,"ctl_water3"],
pch=19,col="black",xlab="ctl_water1",ylab="ctl_water3",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"ctl_water1"],
normalized.gene.expression[,"ctl_water3"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"ctl_water2"],
y = normalized.gene.expression[,"ctl_water3"],
pch=19,col="black",xlab="ctl_water2",ylab="ctl_water3",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"ctl_water2"],
normalized.gene.expression[,"ctl_water3"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"ctl_drought1"],
y = normalized.gene.expression[,"ctl_drought2"],
pch=19,col="black",xlab="ctl_drought1",ylab="ctl_drought2",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"ctl_drought1"],
normalized.gene.expression[,"ctl_drought2"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"ctl_drought1"],
y = normalized.gene.expression[,"ctl_drought3"],
pch=19,col="black",xlab="ctl_drought1",ylab="ctl_drought3",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"ctl_drought1"],
normalized.gene.expression[,"ctl_drought3"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"ctl_drought2"],
y = normalized.gene.expression[,"ctl_drought3"],
pch=19,col="black",xlab="ctl_drought2",ylab="ctl_drought3",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"ctl_drought2"],
normalized.gene.expression[,"ctl_drought3"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"h1_water1"],
y = normalized.gene.expression[,"h1_water2"],
pch=19,col="black",xlab="h1_water1",ylab="h1_water2",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"h1_water1"],
normalized.gene.expression[,"h1_water2"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"h1_water1"],
y = normalized.gene.expression[,"h1_water3"],
pch=19,col="black",xlab="h1_water1",ylab="h1_water3",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"h1_water1"],
normalized.gene.expression[,"h1_water3"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"h1_water2"],
y = normalized.gene.expression[,"h1_water3"],
pch=19,col="black",xlab="h1_water2",ylab="h1_water3",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"h1_water2"],
normalized.gene.expression[,"h1_water3"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"h1_drought1"],
y = normalized.gene.expression[,"h1_drought2"],
pch=19,col="black",xlab="h1_drought1",ylab="h1_drought2",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"h1_drought1"],
normalized.gene.expression[,"h1_drought2"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"h1_drought1"],
y = normalized.gene.expression[,"h1_drought3"],
pch=19,col="black",xlab="h1_drought1",ylab="h1_drought3",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"h1_drought1"],
normalized.gene.expression[,"h1_drought3"]),
digits = 2),
"%"), collapse=""))
plot(x = normalized.gene.expression[,"h1_drought2"],
y = normalized.gene.expression[,"h1_drought3"],
pch=19,col="black",xlab="h1_drought2",ylab="h1_drought3",cex=0.5)
text(x=3,y=12,
labels = paste(c(
"cor = ",
round(100*cor(normalized.gene.expression[,"h1_drought2"],
normalized.gene.expression[,"h1_drought3"]),
digits = 2),
"%"), collapse=""))
En los scatter plots de la similitud de réplicas los resultados han indicado podemos destacar que la correlación entre las réplicas es muy alta, la mayoría alrededor del 99%, indicando así buena relación entre las réplicas.
Luego, podemos ver en la mayoría de gráficos muy baja dispersión entre los puntos formando así en la mayoría de los casos una figura similar a una línea recta.
Por último, también destacar que la réplica ctl_drought2 tiene una mayor dispersión de datos y por lo tanto mayor variabilidad, dado que a la hora de hacer el estudio se observó que ctl_drought2 se distanciaba más con respecto al resto de datos.
Análisis de componentes principales y clustering jerárquico
library(FactoMineR)
library(factoextra)
pca.gene.expression <- data.frame(colnames(normalized.gene.expression),
t(normalized.gene.expression))
colnames(pca.gene.expression)[1] <- "Sample"
res.pca <- PCA(pca.gene.expression, graph = FALSE,scale.unit = TRUE,quali.sup = 1 )
res.hcpc <- HCPC(res.pca, graph=FALSE,nb.clust = 4)
fviz_dend(res.hcpc,k=,
cex = 0.75,
palette = "Set1",
rect = TRUE, rect_fill = TRUE,
rect_border = "Set1",
type="rectangle",
labels_track_height = 1400
)
En el dendograma, en un principio deberíamos tener cada condición agrupada con sus tres réplicas, sin embargo, vemos que el ctl_drought2 está aislado en un cluster porque sus datos se alejan de su grupo condición. Con lo cual podríamos eliminar dicha réplica ya que altera la estabilidad de los datos.
Además, también se visualiza una alteración en el cluster referente a “water”, en el que h1_water3 está desviado de su grupo condición h1_water. Lo mismo ocurre con ctl_water1 que está más cercano al grupo h1_water. Esto es debido a la variabilidad de datos obtenida en estas condiciones.
fviz_pca_ind(res.pca,
col.ind = c("ctl_water", "ctl_water", "ctl_water", "ctl_drought", "ctl_drought", "ctl_drought", "h1_water", "h1_water", "h1_water", "h1_drought", "h1_drought", "h1_drought"),
pointsize = 2,
pointshape = 21,
fill = "black",
repel = TRUE,
addEllipses = TRUE,
ellipse.type = "confidence",
legend.title = "Conditions",
title = "",
show.legend = TRUE,
show.guide = TRUE
)
En el gráfico de elipses de confianza, las elipses indican la agrupación de los datos, en este caso todas las réplicas aparecen agrupadas en sus respectivas condiciones. Tanto ctl_water, como h1_drought como h1_water aparecen correctamente agrupados. Las réplicas de estos grupos no están muy alejadas del punto medio indicando poca variabilidad de los datos.
De nuevo, vemos que la condición ctl_drought2 está bastante dispersa con respecto a las réplicas de su misma condición, formando una elipse que se asemeja más a una recta.
A diferencia del dendograma, vemos en este gráfico que las réplicas están bien agrupadas.
Comparación del transcriptoma generado con el de referencia
Antes de seguir avanzando con el análisis de datos, es conveniente confirmar que la anotación utilizada (es decir, la descargada de la base de datos) es adecuada para el estudio. Para ello, podemos comparar el transcriptoma de la base de datos con los transcritos que se han detectado en total en el estudio (en todas las muestras).
Como podemos ver, el
transcriptoma utilizado contiene el 100% de los transcritos
observados.
Visualización del mapeo
Podemos observar el resultado del mapeo en IGV. En la siguiente imagen se muestran algunas de las lecturas que han mapeado en el gen AT4G33100.
Cada línea representa una lectura, de forma que los rectángulos grises son secuencias del genoma que se encuentran en la lectura y las líneas entre los rectángulos corresponden a fragmentos que están presentes en el genoma, pero no en las lecturas. Como podemos observar en la anotación, estas zonas se corresponden con los intrones del gen (líneas azules).
En primer lugar, vamos a usar scatterplots para visualizar el efecto de los distintos tratamientos que se realizan en el estudio.
# Calculamos la expresión media de cada gen en cada condición
ctl_water <- (normalized.gene.expression[,"ctl_water1"] + normalized.gene.expression[,"ctl_water2"]+ normalized.gene.expression[,"ctl_water3"])/3
ctl_drought <- (normalized.gene.expression[,"ctl_drought1"] + normalized.gene.expression[,"ctl_drought2"]+ normalized.gene.expression[,"ctl_drought3"])/3
h1_water <- (normalized.gene.expression[,"h1_water1"] + normalized.gene.expression[,"h1_water2"]+ normalized.gene.expression[,"h1_water3"])/3
h1_drought <- (normalized.gene.expression[,"h1_drought1"] + normalized.gene.expression[,"h1_drought2"]+ normalized.gene.expression[,"h1_drought3"])/3
mean.expression <- matrix(c(ctl_water,ctl_drought, h1_water,h1_drought),ncol=4)
colnames(mean.expression) <- c("ctl_water","ctl_drought","h1_water","h1_drought")
rownames(mean.expression) <- rownames(normalized.gene.expression)
par(mfrow=c(2,3))
plot(ctl_water,ctl_drought,pch=19,cex=0.7,xlab="ctl_water",
ylab="ctl_drought",cex.lab=1.25,col="black")
plot(h1_water,h1_drought,pch=19,cex=0.7,xlab="h1_water",
ylab="h1_drought",cex.lab=1.25,col="black")
plot(ctl_drought,h1_drought,pch=19,cex=0.7,xlab="ctl_drought",
ylab="h1_drought",cex.lab=1.25,col="black")
plot(ctl_water,h1_water,pch=19,cex=0.7,xlab="ctl_water",
ylab="h1_water",cex.lab=1.25,col="black")
plot(ctl_water,h1_drought,pch=19,cex=0.7,xlab="ctl_water",
ylab="h1_drought",cex.lab=1.25,col="black")
Podemos ver que el transcriptoma no se ve fuertemente alterado por los tratamientos de sequía e inoculación con Aeromonas sp., por lo que tendremos que ser muy permisivos a la hora de determinar qué genes se expresan diferencialmente. En este caso, consideramos genes diferencialmente expresados (DEGs) aquellos genes en los que el cambio observado en la expresión génica sea sustancial y significativo. Se usarán umbrales de |log2fc|>log2(1.5) y q-valor<0.01.
library(limma)
# Definimos el diseño experimental y determinamos las comparaciones que se realizarán
experimental.design <- model.matrix(~ -1+factor(c(1,1,1,2,2,2,3,3,3,4,4,4)))
colnames(experimental.design) <- c("ctl_water","ctl_drought","h1_water","h1_drought")
linear.fit <- lmFit(normalized.gene.expression, experimental.design)
contrast.matrix <- makeContrasts(ctl_drought-ctl_water, h1_drought-h1_water, h1_drought-ctl_drought,h1_water-ctl_water, h1_drought-ctl_water,levels=c("ctl_water","ctl_drought","h1_water","h1_drought"))
contrast.linear.fit <- contrasts.fit(linear.fit, contrast.matrix)
contrast.results <- eBayes(contrast.linear.fit)
Comparación 1. CTL_Drought vs CTL_Water
El objetivo de esta primera comparación es identificar los genes involucrados en la respuesta a sequía.
ctl_drought.ctl_water <- topTable(contrast.results, number=32833,coef=1,sort.by="logFC")
log.fold.change <- ctl_drought.ctl_water$logFC
q.value <- ctl_drought.ctl_water$adj.P.Val
genes.ids <- rownames(ctl_drought.ctl_water)
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c1 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
repressed.c1 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
length(activated.c1)
## [1] 112
length(repressed.c1)
## [1] 36
log.q.val <- -log10(q.value)
plot(log.fold.change,log.q.val,pch=19,col="black",cex=0.8,
xlim=c(-6,6),ylim = c(0,4),
xlab="log2(Fold-chage)",ylab="-log10(q-value)",cex.lab=1.5,
main="CTL_drought vs CTL_water")
points(x = log.fold.change[activated.c1],
y = log.q.val[activated.c1],col="red",cex=0.8,pch=19)
points(x = log.fold.change[repressed.c1],
y = log.q.val[repressed.c1],col="blue",cex=0.8,pch=19)
arrows(log.fold.change["AT5G02020"],log.q.val["AT5G02020"],2.5,1.2,length=0)
text(2.8,1.2,"SIS",cex=0.7)
arrows(log.fold.change["AT4G17500"],log.q.val["AT4G17500"],2,2,length=0)
text(2.5,2.3,"ERF-1",cex=0.7)
arrows(log.fold.change["AT4G08950"],log.q.val["AT4G08950"],-2,2,length=0)
text(-2,2.3,"EX0",cex=0.7)
write.table(activated.c1,file="activated_ctl_drought_ctl_water.txt",row.names = FALSE,quote=FALSE,col.names = FALSE)
write.table(repressed.c1,file="repressed_ctl_drought_ctl_water.txt",row.names = FALSE,quote=FALSE,col.names = FALSE)
Genes activados comparación 1 Genes reprimidos comparación 1
En esta comparación, se han detectado 112 genes activados y 36 genes reprimidos. Entre los genes activados, encontramos SIS (Salt Induced Serine Rich), proteína que se sabe participa en la tolerancia a salinidad (condición estrechamente relacionada a la baja disponibilidad hídrica).
Comparación 2. H1_drought vs H1_water
El objetivo en este caso es determinar qué genes participan en la respuesta a sequía cuando las plantas han sido inoculadas con Aeromonas.
h1_drought.h1_water <- topTable(contrast.results, number=32833,coef=2,sort.by="logFC")
log.fold.change <- h1_drought.h1_water$logFC
q.value <- h1_drought.h1_water$adj.P.Val
genes.ids <- rownames(h1_drought.h1_water)
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c2 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
repressed.c2 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
length(activated.c2)
## [1] 282
length(repressed.c2)
## [1] 45
log.q.val <- -log10(q.value)
plot(log.fold.change,log.q.val,pch=19,col="black",cex=0.8,
xlim=c(-6,6),ylim = c(0,4),
xlab="log2(Fold-chage)",ylab="-log10(q-value)",cex.lab=1.5,
main="H1_drought vs H1_water")
points(x = log.fold.change[activated.c2],
y = log.q.val[activated.c2],col="red",cex=0.8,pch=19)
points(x = log.fold.change[repressed.c2],
y = log.q.val[repressed.c2],col="blue",cex=0.8,pch=19)
arrows(log.fold.change["AT3G23550"],log.q.val["AT3G23550"],2.5,2.5,length=0)
text(2.9,2.5,"DTX18",cex=0.7)
arrows(log.fold.change["AT1G73500"],log.q.val["AT1G73500"],2.3,1.8,length=0)
text(2.6,1.8,"MKK9",cex=0.7)
arrows(log.fold.change["AT5G15600"],log.q.val["AT5G15600"],-2,1,length=0)
text(-2,1,"EX0",cex=0.7)
write.table(activated.c2,file="activated_h1_drought_h1_water.txt",row.names = FALSE,quote=FALSE,col.names = FALSE)
write.table(repressed.c2,file="repressed_h1_drought_h1_water.txt",row.names = FALSE,quote=FALSE,col.names = FALSE)
Genes activados comparación 2 Genes reprimidos comparación 2
Con los criterios seleccionados, encontramos 282 genes activados y 45 genes reprimidos. En este caso, uno de los genes activados es la proteína quinasa activada por mitógenos MKK9, que tiene función de señalización intracelular.
Comparación 3. H1_drought vs CTL_drought.
Con esta comparación, analizamos el efecto que tiene inocular con Aeromonas las raíces de Arabidopsis, en condiciones de sequía. Así, buscamos genes cuya expresión cambia gracias a Aeromonas.
h1_drought.ctl_drought <- topTable(contrast.results, number=32833,coef=3,sort.by="logFC")
log.fold.change <- h1_drought.ctl_drought$logFC
q.value <- h1_drought.ctl_drought$adj.P.Val
genes.ids <- rownames(h1_drought.ctl_drought)
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c3 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
repressed.c3 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
length(activated.c3)
## [1] 84
length(repressed.c3)
## [1] 28
log.q.val <- -log10(q.value)
plot(log.fold.change,log.q.val,pch=19,col="black",cex=0.8,
xlim=c(-6,6),ylim = c(0,4),
xlab="log2(Fold-chage)",ylab="-log10(q-value)",cex.lab=1.5,
main="CTL_drought vs H1_drought")
points(x = log.fold.change[activated.c3],
y = log.q.val[activated.c3],col="red",cex=0.8,pch=19)
points(x = log.fold.change[repressed.c3],
y = log.q.val[repressed.c3],col="blue",cex=0.8,pch=19)
arrows(log.fold.change["AT3G15500"],log.q.val["AT3G15500"],2.5,2.5,length=0)
text(2.9,2.5,"NAC3",cex=0.7)
write.table(activated.c3,file="activated_h1_drought_ctl_drought.txt",row.names = FALSE,quote=FALSE,col.names = FALSE)
write.table(repressed.c3,file="repressed_h1_drought_ctl_drought.txt",row.names = FALSE,quote=FALSE,col.names = FALSE)
Genes activados comparación 3 Genes reprimidos comparación 3
En esta comparación, se detectaron 84 genes activados y 28 genes reprimidos. Entre los genes activados, se incluye NAC3, un factor de transcripción relacionado con procesos de desarrollo.
Comparación 4. H1_water vs CTL_water
Con esta cuarta comparación, estudiamos el efecto de la inoculación con Aeromonas cuando Arabidopsis se encuentra en condiciones basales.
h1_water.ctl_water <- topTable(contrast.results, number=32833,coef=4,sort.by="logFC")
log.fold.change <- h1_water.ctl_water$logFC
q.value <- h1_water.ctl_water$adj.P.Val
genes.ids <- rownames(h1_water.ctl_water)
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c4 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
repressed.c4 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
length(activated.c4)
## [1] 0
length(repressed.c4)
## [1] 0
log.q.val <- -log10(q.value)
plot(log.fold.change,log.q.val,pch=19,col="black",cex=0.8,
xlim=c(-6,6),ylim = c(0,4),
xlab="log2(Fold-chage)",ylab="-log10(q-value)",cex.lab=1.5,
main="H1_water vs CTL_water")
points(x = log.fold.change[activated.c4],
y = log.q.val[activated.c4],col="red",cex=0.8,pch=19)
points(x = log.fold.change[repressed.c4],
y = log.q.val[repressed.c4],col="blue",cex=0.8,pch=19)
Podemos observar que ningún gen se expresa diferencialmente en esta comparación. Cuando las plantas no están sometidas a estrés por sequía, el tratamiento con Aeromonas no tiene ningún efecto sobre el transcriptoma de la parte aérea.
Comparación 5. H1_drought vs CTL_water
Por último, comparamos la expresión génica en las plantas control (regadas y sin inocular) con plantas inoculadas con Aeromonas y en condiciones de sequía. De esta forma, podremos observar un efecto sinérgico entre el tratamiento de sequía y el efecto que produce la bacteria.
h1_drought.ctl_water <- topTable(contrast.results, number=32833,coef=5,sort.by="logFC")
log.fold.change <- h1_drought.ctl_water$logFC
q.value <- h1_drought.ctl_water$adj.P.Val
genes.ids <- rownames(h1_drought.ctl_water)
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c5 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
repressed.c5 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
length(activated.c5)
## [1] 556
length(repressed.c5)
## [1] 150
log.q.val <- -log10(q.value)
plot(log.fold.change,log.q.val,pch=19,col="black",cex=0.8,
xlim=c(-6,6),ylim = c(0,4),
xlab="log2(Fold-chage)",ylab="-log10(q-value)",cex.lab=1.5,
main="H1_drought vs CTL_water")
points(x = log.fold.change[activated.c5],
y = log.q.val[activated.c5],col="red",cex=0.8,pch=19)
points(x = log.fold.change[repressed.c5],
y = log.q.val[repressed.c5],col="blue",cex=0.8,pch=19)
arrows(log.fold.change["AT1G74930"],log.q.val["AT1G74930"],3.8,3.6,length=0)
text(4.3,3.7,"ORA47",cex=0.7)
arrows(log.fold.change["AT5G47220"],log.q.val["AT5G47220"],4,3,length=0)
text(4,2.8,"ERF2",cex=0.7)
write.table(activated.c5,file="activated_h1_drought_ctl_water.txt",row.names = FALSE,quote=FALSE,col.names = FALSE)
write.table(repressed.c5,file="repressed_h1_drought_ctl_water.txt",row.names = FALSE,quote=FALSE,col.names = FALSE)
Genes activados comparación 5 Genes reprimidos comparación 5
En esta última comparación, se han identificado 556 genes activados y 150 genes reprimidos. Se destaca el gen ORA47, factor de transcripción que también pertenece a la familia conocida como ERF/AP2.
Análisis con DESeq2
library(DESeq2)
gene.count.matrix <- read.table(file = "gene_count_matrix.csv",header = T,sep = ",")
gene.ids <- sapply(X = strsplit(x = gene.count.matrix$gene_id,split = "\\|"),FUN = function(x){return(x[1])})
gene.count.matrix <- gene.count.matrix[,-1]
rownames(gene.count.matrix) <- gene.ids
dds <- DESeqDataSetFromMatrix(countData=gene.count.matrix, colData=pheno.data, design = ~ aeromonas + drought)
dds$aeromonas <- relevel(dds$aeromonas, ref = "ctl")
dds$drought <- relevel(dds$drought, ref = "ctl")
dds <- DESeq(dds)
mod_mat <- model.matrix(design(dds), colData(dds))
ctl_water <- colMeans(mod_mat[dds$aeromonas == "ctl" & dds$drought == "ctl", ])
ctl_drought <- colMeans(mod_mat[dds$aeromonas == "ctl" & dds$drought == "drought", ])
h1_water <- colMeans(mod_mat[dds$aeromonas == "h1" & dds$drought == "ctl", ])
h1_drought <- colMeans(mod_mat[dds$aeromonas == "h1" & dds$drought == "drought", ])
# Comparación 1
res1 <- results(dds, contrast = ctl_drought - ctl_water)
genes.ids <- rownames(res1)
log.fold.change <- res1$log2FoldChange
q.value <- res1$padj
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c1.deseq2 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
activated.c1.deseq2 <- activated.c1.deseq2[!is.na(activated.c1.deseq2)]
repressed.c1.deseq2 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
repressed.c1.deseq2 <- repressed.c1.deseq2[!is.na(repressed.c1.deseq2)]
length(activated.c1.deseq2)
## [1] 324
length(repressed.c1.deseq2)
## [1] 183
# Comparación 2
res2 <- results(dds, contrast = h1_drought - h1_water)
genes.ids <- rownames(res2)
log.fold.change <- res2$log2FoldChange
q.value <- res2$padj
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c2.deseq2 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
activated.c2.deseq2 <- activated.c2.deseq2[!is.na(activated.c2.deseq2)]
repressed.c2.deseq2 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
repressed.c2.deseq2 <- repressed.c2.deseq2[!is.na(repressed.c2.deseq2)]
length(activated.c2.deseq2)
## [1] 324
length(repressed.c2.deseq2)
## [1] 183
# Comparación 3
res3 <- results(dds, contrast = h1_drought - ctl_drought)
genes.ids <- rownames(res3)
log.fold.change <- res3$log2FoldChange
q.value <- res3$padj
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c3.deseq2 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
activated.c3.deseq2 <- activated.c3.deseq2[!is.na(activated.c3.deseq2)]
repressed.c3.deseq2 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
repressed.c3.deseq2 <- repressed.c3.deseq2[!is.na(repressed.c3.deseq2)]
length(activated.c3.deseq2)
## [1] 101
length(repressed.c3.deseq2)
## [1] 30
# Comparación 4
res4 <- results(dds, contrast = h1_water - ctl_water)
genes.ids <- rownames(res4)
log.fold.change <- res4$log2FoldChange
q.value <- res4$padj
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c4.deseq2 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
activated.c4.deseq2 <- activated.c4.deseq2[!is.na(activated.c4.deseq2)]
repressed.c4.deseq2 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
repressed.c4.deseq2 <- repressed.c4.deseq2[!is.na(repressed.c4.deseq2)]
length(activated.c4.deseq2)
## [1] 101
length(repressed.c4.deseq2)
## [1] 30
# Comparación 5
res5 <- results(dds, contrast = h1_drought - ctl_water)
genes.ids <- rownames(res5)
log.fold.change <- res5$log2FoldChange
q.value <- res5$padj
names(log.fold.change) <- genes.ids
names(q.value) <- genes.ids
activated.c5.deseq2 <- genes.ids[log.fold.change > log2(1.5) & q.value < 0.1]
activated.c5.deseq2 <- activated.c5.deseq2[!is.na(activated.c5.deseq2)]
repressed.c5.deseq2 <- genes.ids[log.fold.change < - log2(1.5) & q.value < 0.1]
repressed.c5.deseq2 <- repressed.c5.deseq2[!is.na(repressed.c5.deseq2)]
length(activated.c5.deseq2)
## [1] 613
length(repressed.c5.deseq2)
## [1] 282
Comparación resultados limma vs DESeq2
library(ggvenn)
ggvenn(list(limma=c(activated.c2,repressed.c2),DESeq2=c(activated.c2.deseq2,repressed.c2.deseq2)))
ggvenn(list(limma=c(activated.c5,repressed.c5),DESeq2=c(activated.c5.deseq2,repressed.c5.deseq2)))
En estos diagramas de Venn, comparamos los DEGs obtenidos al usar limma con los identificados por DESeq2 en las comparaciones h1_drought vs h1_water y h1_drought vs ctl_drought.
Como podemos ver, los conjuntos de genes obtenidos al analizar los datos con limma y con DESeq2 son bastante diferentes. No esperamos un solapamiento total, ya que existen diferencias entre los análisis con limma (que usa fpkm para medir la expresión génica) y DESeq2 (conteos). Como es típico, DESeq2 identifica un mayor número de DEGs, al realizar análisis más sensibles. Sin embargo, DESeq2 también suele dar lugar a un mayor número de falsos positivos.
Además, en este estudio en concreto, hemos observado que el efecto sobre el transcriptoma es relativamente sutil. En caso de tener un efecto más acusado, en el que pudiéramos establecer umbrales de fold-change y de q-valores más restrictivos, puede que observáramos una concordancia mejor entre los análisis.
Comparaciones entre conjuntos de genes
ggvenn(list(degs_drought=c(activated.c1,repressed.c1),degs_h1_drought=c(activated.c5,repressed.c5)))
En este diagrama de Venn, comparamos dos conjuntos de genes:
DEGs Ctl_water vs Ctl_drought
DEGs Ctl_water vs H1_drought
Podemos observar que existen 116 genes que cambian en ambos casos, es decir, genes que son activados por las condiciones de sequía pero en los que no está implicado Aeromonas. Por otro lado, existen 590 genes que cambian al aplicar la sequía a plantas inoculadas pero no en plantas sin inocular. Estos genes deben de ser aquellos cuya expresión cambia gracias a Aeromonas, y confieren a la planta una mejor tolerancia a la sequía.
La ontología de genes permite la incorporación de información de forma sistemática e inequívoca a genes mediante anotación. Esta se compone de un vocabulario estructurado (de más genérico a más específico) y controlado de términos que describen los productos génicos en términos de procesos biológicos, componentes celulares y funciones moleculares.
Un enriquecimiento de términos de ontología de genes tiene como objetivo mejorar la representación de la información genética, haciéndola más precisa y completa, de forma que sea posible interpretar de manera efectiva los resultados.
library(clusterProfiler)
library(org.At.tair.db)
library(enrichplot)
Genes diferencialmente activados H1_drought vs CTL_drought
En este caso se realiza un enriquecimiento de términos GO del conjunto de genes diferencialemente activados en las plantas de Arabidopsis thaliana inoculadas con Aeromonas sp. en condiciones de sequía en comparación con las plantas no inoculadas en condiciones de sequía.
activated_h1_d_ctl_d <- read.table(file="activated_h1_drought_ctl_drought.txt")[[1]]
activated_h1_d_ctl_d.atha.enrich.go <- enrichGO(gene = activated_h1_d_ctl_d,
OrgDb = org.At.tair.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = FALSE,
keyType = "TAIR")
df.1 <- as.data.frame(activated_h1_d_ctl_d.atha.enrich.go)
write.table(as.data.frame(df.1$pvalue,df.1$ID), file = "activated_h1_d_ctl_d_revigo", quote = FALSE)
| Término de GO y q-valor | Descripción | Genes Representativos |
|---|---|---|
| GO:0002376 7.369807e-05 | immune system process | ERF2/PEN3/CYP81F2 |
| GO:0008219 0.0005994776 | cell death | PEN3/PLA2A/NUDT7 |
| GO:0009867 1.141716e-12 | jasmonic acid mediated signaling pathway | JAZ7/ERF2/NAC3 |
| GO:0042430 3.081945e-05 | indole-containing compound metabolic process | HSPRO2/PEN3/CYP81F2 |
| GO:0042908 0.004019064 | xenobiotic transport | DTX18 |
| GO:0090693 0.001098587 | plant organ senescence | BGLU11/MKK9/NAC6 |
En la tabla se puede observar el proceso biológico mayormente representado de cada grupo de procesos relacionados entre sí que se encuentran significativamente alterados. Entre los más representados, el de la ruta de señalización mediada por el ácido jasmónico (JA) está significativamente enriquecido en comparación con los demás al tener el qvalor más bajo. Estas observaciones se pueden corroborar al ver la imagen del treemap en la que se aprecia como el bloque más grande de todos es el correspondiente a esta ruta. Además, el tamaño del mismo grupo (bloques del mismo color) ocupa claramente casi toda la imagen indicando que son procesos muy importantes para la resistencia al estrés hídrico proporcionada por Aeromonas sp., es por eso que son rutas mayormente activadas.
Para más inri, el análisis GO destaca DEGs relacionados con JA, como JAZ7, ERF2 y NAC3, que se sabe que promueven la resistencia de las plantas al estrés por deshidratación. NAC3 se ha obtenido anteriormente en los análisis de genes diferencialmente activados para estas condiciones, y ERF2 parece estar también involucrado en otros procesos enriquecidos. También podemos destacar PEN3 al regular varios proceso de importancia.
Genes diferencialmente reprimidos H1_drought vs CTL_drought
En este caso se realiza un enriquecimiento de términos GO el conjunto de genes diferencialemente reprimidos en las plantas de Arabidopsis thaliana inoculadas con Aeromonas sp. en condiciones de sequía en comparación con las plantas no inoculadas en condiciones de sequía.
repressed_h1_d_ctl_d <- read.table(file="repressed_h1_drought_ctl_drought.txt")[[1]]
repressed_h1_d_ctl_d.atha.enrich.go <- enrichGO(gene = repressed_h1_d_ctl_d,
OrgDb = org.At.tair.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = FALSE,
keyType = "TAIR")
df.2 <- as.data.frame(repressed_h1_d_ctl_d.atha.enrich.go)
write.table(as.data.frame(df.2$pvalue,df.2$ID), file="repressed_h1_d_ctl_d_revigo", quote=FALSE)
| Término de GO y q-valor | Descripción | Genes Representativos |
|---|---|---|
| GO:0009639 0.0002585459 | response to red or far red light | ERF34/ERD9/BBX32 |
| GO:0009813 0.009531447 | flavonoid biosynthetic process | TT4/FLS1 |
| GO:0031537 0.0118085 | regulation of anthocyanin metabolic process | TT4 |
En este caso es posible ver tanto en el treemap como en la tabla que hay un enriquecimiento de genes reprimidos en respuesta a la luz roja o roja lejana en las plantas inoculadas, seguido del proceso de biosíntesis de flavonoides.
Genes diferencialmente activados H1_drought vs CTL_water
En este caso se realiza un enriquecimiento de términos GO el conjunto de genes diferencialmente activados en las plantas de Arabidopsis thaliana inoculadas con Aeromonas sp. en condiciones de sequía en comparación con las plantas no inoculadas y regadas.
activated_h1_d_ctl_w <- read.table(file="activated_h1_drought_ctl_water.txt")[[1]]
activated_h1_d_ctl_w.atha.enrich.go <- enrichGO(gene = activated_h1_d_ctl_w,
OrgDb = org.At.tair.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = FALSE,
keyType = "TAIR")
df.3 <- as.data.frame(activated_h1_d_ctl_w.atha.enrich.go)
write.table(as.data.frame(df.3$pvalue,df.3$ID), file="activated_h1_d_ctl_w_revigo", quote=FALSE)
| Término de GO y q-valor | Descripción | Genes Representativos |
|---|---|---|
| GO:0002376 4.074928e-18 | immune system process | ERF2/PLA2A/PEN3 |
| GO:0008219 1.272887e-08 | cell death | PLA2A/PEN3/NAC6 |
| GO:0009751 2.481618e-49 | response to salicylic acid | BT2/ATS40.4/SYP122 |
| GO:0009966 1.864563e-14 | regulation of signal transduction | JAZ7/JAZ1/JAZ6 |
| GO:0010150 3.679929e-21 | leaf senescence | RD26/MKK9/SAG21 |
| GO:0012501 8.360128e-09 | programmed cell death | PLA2A/PEN3/NAC6 |
| GO:0034219 0.003716663 | carbohydrate transmembrane transport | ERD6/STP1/ESL1 |
| GO:0042430 2.39423e-23 | indole-containing compound metabolic process | IGMT1/UGT74E2/MKK9 |
| GO:0045229 0.01758479 | external encapsulating structure organization | PEN3/TCH4/CYP81F2 |
Cabe resaltar la respuesta al ácido salicílico por el gran enriquecimiento que tiene en comparación con el resto de procesos, al igual que todo el grupo en comparación con el resto de grupos.
En relación a los genes más representativos se puede destacar el ERF2 correspondiente a la activación del sistema inmunitario, uno de los procesos mayormente enriquecidos. Por otro lado, volvemos a ver el gen JAZ7 al igual que en la comparación de genes activados H1_drought vs CTL_drought, lo que da a entender que es fundamental para la respuesta al estrés hídrico.
Genes diferencialmente reprimidos H1_drought vs CTL_water
En este caso se realiza un enriquecimiento de términos GO el conjunto de genes diferencialmente reprimidos en las plantas de Arabidopsis thaliana inoculadas con Aeromonas sp. en condiciones de sequía en comparación con las plantas no inoculadas y regadas.
repressed_h1_d_ctl_w <- read.table(file="repressed_h1_drought_ctl_water.txt")[[1]]
repressed_h1_d_ctl_w.atha.enrich.go <- enrichGO(gene = repressed_h1_d_ctl_w,
OrgDb = org.At.tair.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = FALSE,
keyType = "TAIR")
df.4 <- as.data.frame(repressed_h1_d_ctl_w.atha.enrich.go)
write.table(as.data.frame(df.4$pvalue,df.4$ID), file="repressed_h1_d_ctl_w_revigo", quote=FALSE)
| Término de GO y q-valor | Descripción | Genes Representativos |
|---|---|---|
| GO:0000278 0.001719728 | mitotic cell cycle | CYCA1;1/TPX2/CYCB1;1 |
| GO:0007017 0.01155884 | microtubule-based process | TPX2/SP1L4/PAKRP2 |
| GO:0008356 0.01895755 | asymmetric cell division | TMM/BASL |
| GO:0009411 0.0008523922 | response to UV | ELIP1/MYB4/F3H |
| GO:0009813 0.0001719778 | flavonoid biosynthetic process | FLS1/F3H/RHM1 |
| GO:0010374 0.001380033 | stomatal complex development | TMM/CKX6/EPF2 |
| GO:0042546 0.004344226 | cell wall biogenesis | ERF34/AGP17/RGP1 |
| GO:0051338 0.003546343 | regulation of transferase activity | CYCA1;1/TPX2/CDC20.1 |
En este caso no se observan grandes diferencias en cuanto al enriquecimiento de procesos o grupos de procesos siendo los q-valores bastante parecidos y el treemap bastante homogéneo en cuanto a tamaños. Aunque los más enriquecidos coinciden con los de la condición de los genes diferencialmente reprimidos H1_drought vs CTL_drought.
Enriquecimiento de términos de rutas KEGG
# Genes diferencialmente activados H1_drought vs CTL_drought
activated_h1_d_ctl_d.atha.enrich.kegg <- enrichKEGG(gene = activated_h1_d_ctl_d,
organism = "ath",
pAdjustMethod = "BH",
pvalueCutoff = 0.05)
df.activated_h1_d_ctl_d.atha.enrich.kegg <- as.data.frame(activated_h1_d_ctl_d.atha.enrich.kegg)
df.activated_h1_d_ctl_d.atha.enrich.kegg[,1:4]
## category subcategory ID
## ath04075 Environmental Information Processing Signal transduction ath04075
## ath00592 Metabolism Lipid metabolism ath00592
## Description
## ath04075 Plant hormone signal transduction - Arabidopsis thaliana (thale cress)
## ath00592 alpha-Linolenic acid metabolism - Arabidopsis thaliana (thale cress)
# Genes diferencialmente activados H1_drought vs CTL_water
activated_h1_d_ctl_w.atha.enrich.kegg <- enrichKEGG(gene = activated_h1_d_ctl_w,
organism = "ath",
pAdjustMethod = "BH",
pvalueCutoff = 0.05)
df.activated_h1_d_ctl_w.atha.enrich.kegg <- as.data.frame(activated_h1_d_ctl_w.atha.enrich.kegg)
df.activated_h1_d_ctl_w.atha.enrich.kegg[,1:4]
## category subcategory ID
## ath00592 Metabolism Lipid metabolism ath00592
## ath04626 Organismal Systems Environmental adaptation ath04626
## Description
## ath00592 alpha-Linolenic acid metabolism - Arabidopsis thaliana (thale cress)
## ath04626 Plant-pathogen interaction - Arabidopsis thaliana (thale cress)
En ambos conjuntos de genes, existe un enriquecimiento en genes de la ruta del metabolismo del ácido alfa-linolénico. A continuación, visualizaremos esta ruta.
library(pathview)
log.fold.change <- h1_drought.ctl_water$logFC
names(log.fold.change) <- rownames(h1_drought.ctl_water)
log.fold.change <- log.fold.change[c(activated.c5,repressed.c5)]
pathview(gene.data = log.fold.change,
pathway.id = "ath00592",
limit = list(gene=max(abs(log.fold.change)), cpd=1),
low = "red", high="green",
species = "ath", gene.idtype = "TAIR")
## [1] "Note: 14 of 706 unique input IDs unmapped."
En esta representación, los genes coloreados representan DEGs (en la comparación h1_drought vs ctl_water). Colores más rojizos indican que la expresión génica de esa enzima es menor en h1_drought (log fold change negativos) mientras que los rectángulos verdes representan una activación de la expresión (log fold change positivos) en el tratamiento. En este caso, no existe una inhibición general de la ruta, ni una activación, sino que unas enzimas son activadas y otras reprimidas, de forma que se acumulan determinados metabolitos mientras disminuye la concentración de otros.
# Genes diferencialmente reprimidos H1_drought vs CTL_drought
repressed_h1_d_ctl_w.atha.enrich.kegg <- enrichKEGG(gene = repressed_h1_d_ctl_w,
organism = "ath",
pAdjustMethod = "BH",
pvalueCutoff = 0.05)
df.repressed_h1_d_ctl_w.atha.enrich.kegg <- as.data.frame(repressed_h1_d_ctl_w.atha.enrich.kegg)
df.repressed_h1_d_ctl_w.atha.enrich.kegg[,1:4]
## category subcategory ID
## ath00941 Metabolism Biosynthesis of other secondary metabolites ath00941
## Description
## ath00941 Flavonoid biosynthesis - Arabidopsis thaliana (thale cress)
# Genes diferencialmente reprimidos H1_drought vs CTL_water
repressed_h1_d_ctl_d.atha.enrich.kegg <- enrichKEGG(gene = repressed_h1_d_ctl_d,
organism = "ath",
pAdjustMethod = "BH",
pvalueCutoff = 0.05)
df.repressed_h1_d_ctl_d.atha.enrich.kegg <- as.data.frame(repressed_h1_d_ctl_d.atha.enrich.kegg)
df.repressed_h1_d_ctl_d.atha.enrich.kegg[,1:4]
## category subcategory ID
## ath00941 Metabolism Biosynthesis of other secondary metabolites ath00941
## Description
## ath00941 Flavonoid biosynthesis - Arabidopsis thaliana (thale cress)
En cuanto a los genes reprimidos en condiciones de sequía en plantas inoculadas con Aeromonas, encontramos un enriquecimiento en la ruta de síntesis de los flavonoides.
pathview(gene.data = log.fold.change,
pathway.id = "ath00941",
limit = list(gene=max(abs(log.fold.change)), cpd=1),
low = "red", high="green",
species = "ath", gene.idtype = "TAIR")
## [1] "Note: 14 of 706 unique input IDs unmapped."
Se muestra la ruta de síntesis de los flavonoides, indicando con rectángulos rojizos aquellas enzimas cuya expresión disminuye al aplicar condiciones de sequía e inocular con Aeromonas. En este caso, todos los DEGs identificados corresponden a genes reprimidos, lo que parece indicar que se produce una inhibición general de la síntesis de flavonoides.
Con los análisis realizados, hemos obtenido los genes diferencialmente expresados para condiciones de estrés por sequía y al ser inoculada con Aeromonas sp. Hemos observado cómo Arabidopsis thaliana en condiciones de sequía es capaz de producir una respuesta más fuerte a este estrés abiótico (expresa más genes de forma diferencial) cuando está tratada con Aeromonas sp.
Al estudiar el enriquecimiento de términos GO, hemos podido comprobar que, para nuestra sorpresa, los genes que se activan al inocular con Aeromonas las plantas en condiciones de sequía están involucradas en la respuesta a estrés hídrico. Nuestros análisis de enriquecimiento también han revelado que el ácido jasmónico juega un papel importante en este proceso. Por otro lado, hemos visto cómo al comparar condiciones en las que planta está regada, el tratamiento con Aeromonas sp. no tiene un efecto perjudicial para la misma.
Al igual que en nuestro análisis, al realizar los diagramas de Venn, en el artículo (He et al., 2022) se puede apreciar un comportamiento sinérgico entre Aeromonas sp. y la respuesta a estrés por sequía de Arabidopsis thaliana. Además, en el artículo también apuntan a la importancia de la ruta de señalización mediada por ácido jasmónico y aparecen ciertos genes de especial importancia para la resistencia a la sequía como son JAZ7, ERF2 y NAC3. Por otro lado, en el estudio muestran como los flavonoides tienen un papel importante en la interacción planta-microorganismo. Esto es coherente con el hecho de que los genes de la ruta de síntesis de flavonoides sean reprimidos cuando la planta es inoculada con Aeromonas sp. (como hemos visto en el análisis de enriquecimiento de términos de rutas KEGG).
Un experimento sencillo pero eficaz que nos ayudaría a explicar los efectos que Aeromonas sp tendría sobre Arabidopsis thaliana, sería comparar dos plantas en condiciones de deshidratación, siendo una de ellas tratadas con la bacteria y la otra no. Para ello, transplantaríamos las dos plantas con una semana de vida y cultivadas previamente bajo las mismas condiciones, a un medio exactamente igual. Tras este paso, a una de las plantas de le inocula una cierta cantidad controlada de Aeromonas sp y desde ese mismo momento, comienzan las condiciones de deshidratación. Aproximadamente tras 2 semanas de cultivo, se vuelven a regar. De esta forma se puede comprobar que la planta a la que se le inocula la bacteria, muestra un proceso de recuperación que la planta sin Aeromonas, no consigue alcanzar.
He, D., Singh, S. K., Peng, L., Kaushal, R., Vílchez, J. I., Shao, C., ... & Zhang, H. (2022). Flavonoid-attracted Aeromonas sp. from the Arabidopsis root microbiome enhances plant dehydration resistance. The ISME Journal, 16(11), 2622-2632.